Spatial autocorrelation signatures of ecological determinants on plant community characteristics in high Andean wetlands

Understanding how biological communities are shaped is a central tenet of community ecology. Recent evidence highlights the potential of decoupling diversity spatial autocorrelation into its positive and negative components to reveal community assembly processes that would otherwise remain undetected, as well as to improve understanding of their impacts on different facets of diversity. Yet, such approaches have only been implemented to investigate the effects of a few assembly drivers on a small number of diversity components. Here, we used high Andean wetland plant communities over a strong latitudinal gradient to investigate the effects of various ecological factors on spatial autocorrelation patterns of nine community metrics with different informative values, including measures of richness, dominance, evenness and beta-diversity. By combining Moran’s Eigenvector Maps, partial least squares structural equation modeling, and regression analyses, we revealed two groups of community parameters presenting contrasting spatial patterns due to specific sensitivities to ecological factors. While environmental variation and wetland connectivity increased positive spatial autocorrelation in richness and dominance-related parameters, species co-occurrence promoted negative spatial autocorrelation in evenness-related parameters. These results offer new insights regarding both how ecological processes affect species assembly, as well as the information captured by classical taxonomic parameters.

www.nature.com/scientificreports/ negative components (S − (x)). It has the potential to reveal impacts that would remain undetected without such partitioning. For instance, it revealed that environmental disturbance enhances negative spatial autocorrelation in community diversity 14,15,18 without necessarily affecting overall spatial autocorrelation 15,18 . In addition, it showed that disturbance provokes stronger responses in evenness than in richness 14 . These examples demonstrate that the autocorrelation decomposition framework of Biswas et al. 14 can both offer novel insights into how assembly processes affect community assembly and shed light on the ecological signals that different diversity indices are able to capture. Yet, it has so far only been applied to a few community characteristics and assembly drivers.
Here, we used the autocorrelation decomposition approach to characterize the impacts of a range of ecological factors on spatial patterns of taxonomic community parameters and identify which they alter the most. To do so, we analyzed plant communities from high Andean wetlands in the Norte Chico region in Chile, which offer an excellent model for the purpose of our study. Since these ecosystems are insular formations distributed across a marked aridity gradient, various processes are expected to leave imprints at a range of spatial scales. These include dispersal restriction and environmental filtering, two potentially important assembly mechanisms with predicted broad-scale imprints (Table 1), and ecological drift and plant-plant interactions with projected influences at the local scale (i.e., individual wetlands, Table 1). The effects of dispersal on S + (x) is suggested by the patchy and highly isolated nature of these habitats, which likely allows migration only between wetlands in relatively close geographic proximity, and may thus precipitate stronger similarities in community composition 19 and diversity among closest sites 20 . The broad-scale impacts of environmental filtering are anticipated due to the pronounced latitudinal climate gradient in the region, which gradually transitions from Mediterranean in the south to arid in the north 21,22 . The importance of drift, the random changes in species relative abundances 8 , is substantiated by theory and empirical evidence suggesting that it is a major driver of diversity in island-type habitats 20,23,24 . Its actual importance may depend on landscape factors influencing population size within communities and thereby risks of demographic stochasticity 25 . As for plant-plant interactions, they are well established key drivers of alpine vegetation 26,27 . At high elevations, where assemblages are subjected to harsh climatic conditions, positive interactions in particular are expected to assume high importance 28 . They are promoted by facilitator organisms, such as cushion plants, that act as foundation species 29,30 , stabilizing and improving local abiotic Table 1. Expected spatial signatures of community assembly processes on different community characteristics. S + (x): positive spatial autocorrelation, S − (x): negative spatial autocorrelation.

Assembly process
Expected spatial signature Effects on community characteristics Most prominent expected pattern Dispersal Driver of S + (x). Dispersal is a contagious process, generating similarities among nearby communities Dispersal is expected to increase alphadiversity in general. The theory of Island Biogeography 54 postulates a major effect of dispersal on richness. Analyses of plant, invertebrate and macroinvertebrate communities suggest a direct effect of migration on species richness 44 Dispersal is also a main driver of community assembly increasing similarity among connected communities. Empirical studies indicate a major role of dispersal on community uniqueness [92][93][94] Driver of S + (x) in richness and local contribution to beta-diversity

Environmental filtering
Driver of S + (x) essentially. Environmental filtering can provoke similarities among nearby communities through induced spatial dependence due to the non-random distribution of the environmental variables acting as filters. Such effects are expected to be important along environmental gradients Occasional driver of S − (x) When sharp environmental variations occur (i.e., as expected in high mountains), environmental filtering can drive dissimilarities among nearby communities According to empirical evidence, environmental filtering more strongly influences species richness than evenness 14,41 Environmental heterogeneity is also a common driver of community divergence [95][96][97][98] . A major role of environmental variation on community uniqueness is commonly reported [99][100][101] Driver of S + (x) in richness and LCBD www.nature.com/scientificreports/ conditions (e.g. soil nutrient content), supporting the growth of other plants, affecting species recruitment, and enhancing alpine plant diversity [29][30][31][32] .
To investigate the spatial footprints of these processes on high Andean plant diversity, we analyzed how associated ecological factors (wetland connectivity for dispersal 33 , environmental variation for environmental filtering, wetland size for drift, and species co-occurrence for species interactions 23,25,34 influence spatial patterns in community characteristics. Since assembly processes can alter different aspects of the communities (Table 1), we explored ecological factor effects on various alpha-diversity measures (see Supplementary Table S1 online), including measures of richness, dominance (Shannon entropy, Shannon diversity and Simpson diversity), evenness (Pielou's evenness, Shannon evenness and Simpson evenness), as well as on a beta-diversity parameter (i.e., local contribution to beta-diversity, LCBD 35 ). We postulate that ecological factors associated with broad-scale processes (wetland connectivity and environmental variation) 36 will drive S + (x) 13,16 and have stronger effects on richness-related parameters 37,38 and LCBD (Table 1). In addition, we expect ecological factors associated with fine-scale processes (wetland size and species co-occurrence) to generate S − (x) in LCBD, evenness and dominance-related parameters 14,15 (Table 1). To test these predictions, we analyzed positive and negative autocorrelation patterns of each of the community metrics using Moran Eigenvector Maps. Then, we used partial least squares structural equation modeling (PLS-SEM) to estimate the influence of the ecological factors (abiotic environmental variation, species co-occurrence, wetland connectivity, and size) on each plant community parameters. Finally, we investigated whether the influence of ecological factors drove negative or positive spatial autocorrelation in community parameters by investigating the effects of the PLS-SEM path coefficients after excluding potential confounding effects due to the relatedness of the community parameters.

Results
Spatial structure of the diversity parameters. All diversity indices displayed higher levels of S + (x) than S − (x) (Fig. 1). Substantial variation was revealed for both autocorrelation components, however, as indicated by coefficients of variation (CVs) of 50% and 38.5% for S + and S − (x), respectively. Two groups of alphadiversity indices stood out based on their spatial patterns. The first, comprising N0, N1, H and N2, displayed the highest and only significant levels of S + (x), but were characterized by the lowest levels of S − (x), while the second, comprising E10, E20, J and TB, demonstrated the opposite. Autocorrelation patterns of LCBD were intermediate between the two groups of alpha-diversity metrics. Overall, the community metrics displaying the highest S + (x) tended to be associated with low levels of S − (x), and vice versa, as indicated by the strong negative relationship between the two autocorrelation components (Pearson correlation: r = − 0.90, likelihood ratio test of a mixed model of S + (x) in relation to S − (x) including the correlation matrix of community metrics as a random effect: LR = 16.62, df = 1, P < 0.001). Neither the levels of S + (x) nor those of S − (x) were influenced by the covariance of the community metrics (see Supplementary Table S2   www.nature.com/scientificreports/

Influence of local ecological factors on high Andean plant community characteristics.
Among the four ecological factors, only wetland connectivity and size had consistent effects on diversity indices (Table 2). Connectivity increased all community metrics, and significantly so in the case of N1 and N0. In contrast, the estimated effects of wetland size were negative, indicating that increases in wetland size tended to produce a decrease in the community parameters and significantly so in N0. The impacts of environmental variation and species co-occurrence were inconsistent, with path coefficients ranging from negative to positive and characterized by very large CVs (i.e., 243 and 833% respectively, Table 2). Environmental effects were negative except for N0, N1 and N2 ( Table 2, Fig. 2b). Species co-occurrence effects were positive for richness and dominance-influenced parameters (i.e., H and N1), but negative for evenness metrics (E10 and E20). Overall, wetland connectivity and environmental variation had the strongest effects on the community parameters with mean absolute effects being at least three times higher than those of species co-occurrence and wetland size (e.g., mean of path coefficient absolute values in Table 2).
Effects of local ecological factors on positive and negative spatial autocorrelation components of diversity indices. Three of the four ecological factors increased spatial autocorrelation. Significant relationships between the ecological factor effects (i.e., PLS-SEM path coefficients) and levels of autocorrelation of the community parameters were actually revealed for wetland connectivity, environmental variation and species co-occurrence (Table 3). S + (x) levels showed a significant relationship with environmental variation and wetland connectivity effects, and the S − (x) levels with effects of species co-occurrence. Wetland connectivity and environmental variation both increased S + (x) of the community indices that they positively influenced (Table 3, Fig. 2a,b), and particularly N0, N1, N2 and H (Fig. 2a,b,d). S − (x) levels increased in response to negative effects of species occurrence on community indices, leading to increased levels of S -(x) in E10, E20 and J (Fig. 2a,c).

Discussion
By investigating the effects of ecological determinants on the spatial autocorrelation components of different diversity indices, we were able to reveal spatial footprints of assembly drivers and identify groups of diversity metrics presenting specific responses. These results thus provide new insights both on how ecological processes affect species assembly as well as the information captured by classical taxonomic diversity indices.
Here, we identified a range of drivers of S + (x) and S − (x) in high Andean plant communities. Environmental variation and wetland connectivity were two factors anticipated to promote similarity in community diversity among wetlands in close geographic proximity 19,39 . The fact that these factors had the strongest detected effects on community parameters is consistent with previous evidence that suggests they each play a key role in spatial patterns of high Andean wetland communities and diversity 19,20,39,40 . Their impacts on S + (x) in community metrics is also in line with literature associating dispersal restriction and environmental filtering with broad spatial scales 14,41 when dispersal is a limiting factor in structuring communities 42 and in the presence of environmental gradients 16 , as was the case here. Our results demonstrate that not all diversity components were equally influenced by these factors, however. Environmental variation and wetland connectivity increased similarity among wetlands in close geographic proximity in species richness and dominance-related indices (N0, N2, N1 and H), but not in evenness. This adds to the mounting body of evidence showing that richness and evenness respond differently to a wide range of factors and mechanisms 14,37,41,43-50 . As predicted by theoretical literature and earlier empirical results 14,15,17,41 , we found that habitat connectivity increased species richness and S + (x) levels in Table 2. Influence of wetland connectivity and size, species co-occurrence and environmental variability on each plant community parameter estimated by partial least squares structural equation modeling. The confidence intervals and significance of path coefficients of each ecological factor were evaluated based on 10,000 bootstrapping iterations. Significant path coefficients are in bold. Models were carried out separately for each community metric. N0: Richness; H: Shannon entropy; N1: Shannon diversity; N2: Simpson diversity; E10: Shannon evenness; E20: Simpson evenness; J: Pielou's evenness, TB: Total biomass and LCBD: Local contribution to beta-diversity. www.nature.com/scientificreports/ this diversity component. High levels of connectivity also reduced community dominance (i.e., in N1), which is consistent with theoretical expectations that dominance increases with low migration rates, and vice versa 51 .
A noteworthy finding was that environmental effects promoted different spatial patterns on the community parameters depending on their specific impacts. They increased S + (x) in the community metrics that they positively influenced (i.e., richness and dominance-related parameters), which corroborates our expectations that environmental filtering acts at large spatial scale in our study system. However, and contrary to this expectation, no spatial signatures of environmental effect were revealed on biomass and LCBD, while both were significantly decreased by some abiotic variation. This highlights the existence of non-spatially structured physical variations of importance for the development and composition of high Andean wetland plant communities. High Andean landscapes are highly topographically heterogeneous displaying abrupt changes in local edaphic and microclimatic conditions 52,53 , which may lack clear spatial patterns at the resolution and scale of our study. Slope www.nature.com/scientificreports/ and mean annual wind speed, both ecological variables used in this study, are often associated with large local fluctuations and could thus alter plant community characteristics without leaving a specific spatial footprint. Species co-occurrence was the only factor promoting S − (x) in plant community characteristics. Contrary to our expectations, wetland size had no such effect. The negative relationship between wetland size and richness in fact strongly differs from the premise that larger habitats hold a larger number of species as a consequence of higher levels of habitat diversity and propensity to support species from the regional pool 54,55 . Since our results contrast strongly with these expectations, it is suggested that vegetation extent, used here to assess wetland size, is not an adequate proxy for wetland capacity to support different species. The negative effect on richness instead suggests that it reflects the development of competitively superior species. Regarding the effects of species cooccurrence on S -(x), they were found to promote dissimilarities among geographically close wetlands in evenness indices (E20, J and E10). While species co-occurrences do not provide unequivocal proof of biotic interactions 56 , the detected effects are consistent with empirical evidence that species interactions affect evenness 14,37,44,57 , generate fine-scale patterns in diversity 15,58 and play an important role in alpine plant community assembly 29,59 .
A striking finding of our study is the clear separation of two groups of alpha-diversity parameters presenting contrasting spatial patterns due to their specific sensitivities to ecological factors. Since our analyses discarded any confounding effects of the non-independency of the indices, this result cannot be considered as a simple mathematical artefact of their relationships, and thus demonstrates that assembly mechanisms tend to influence specific aspects of diversity. Understanding the informative value of diversity metrics is a longstanding issue in ecology [60][61][62][63][64][65][66][67] , and debate is ongoing regarding which measures are optimal 65,67-69 , as well as whether a single index could be universally applied for comparative purposes 69,70 . The fact that two groups of diversity metrics bore different spatial footprints demonstrates that each group is sensitive to different assembly processes, and encourages their combined use to aid in a better and more thorough understanding of the ecological processes that sustain community assembly. In addition, it shows that each diversity index holds specific potential for conservation and restoration purposes by motivating for either local or landscape-scale actions. While emerging trends from this and previous studies suggest that evenness tends to be particularly sensitive to local/fine-scale processes and richness to broader scale mechanisms, our understanding of how assembly factors generate and alter spatial patterns of different diversity facets is still limited to a few studies 14,15,18,71 . More empirical evidence from a wider range of ecosystems, including non-patchy habitats, is therefore required to establish the generality of these patterns. This will be essential to determining the specific potential of diversity parameters to decipher processes shaping communities.

Methods
Data collection and production. Sampling and processing of high Andean wetland plant communities. We used plant data of Bertin et al. 20 and Pfeiffer et al. 40 collected in 21 high Andean wetlands along a 600 km latitudinal range (between 26°S and 32°S) in north-central Chile between altitudes of 2852 and 4307 m a.s.l. (See Fig. 2 in Bertin et al. 20 ). High Andean wetlands are permanent features formed by groundwater upwelling and occur in the low Alpine and sub-Alpine belts of the Central Andes 72 , fed by glacial streams, snowmelt, and rainfall 72 . These insular features contrast markedly from the arid grassland matrix in which they are found 73 . In the Norte Chico region of Chile, the location of this study, they occur along a latitudinal climatic gradient that varies from hyper-arid in the north to Mediterranean in the south, with mean annual precipitation extremes ranging from 35 to 200 mm, respectively. They support a wide variety plants and animals, including a number of endemic and rare species 72,74 .
Plant assemblages were characterized for each wetland from five 30 × 30 cm vegetation quadrats 20 . The biomass (g/m 2 ) of each species was estimated by summing dried biomass of all individuals present within a quadrat after complete drying of the plant samples at 70 °C. A total of 52 species belonging to 21 families were identified. The plant database can be found as Supplementary Table S3 online.
Community parameters. Nine community metrics (see Supplementary Table S1 online) were estimated after pooling data from the five vegetation quadrats. They included total biomass estimated as biomass (dry weight in gr/m2) of all species present in the community samples, richness, three dominance-related metrics (Shannon entropy, Shannon diversity; Simpson diversity) and three evenness-related metrics (Shannon evenness, Simpson Table 3. Effects of ecological factors on autocorrelation levels of plant community parameters according to an exhaustive model selection approach based on the Akaike information criterion corrected for small sample size. Ecological factor effects were tested separately for the positive (S + (x)) and negative (S − (x)) spatial autocorrelation components of the community parameters. The models used the autocorrelation levels of each community parameter as the dependent variable and the PLS-SEM path coefficients measuring the direct effects of the ecological factors on each community characteristic as predictors. www.nature.com/scientificreports/ evenness, Pielou's evenness) and the local contribution to beta-diversity calculated from the Hellinger transformed community data using the function beta.div in the adespatial package 35 .
Ecological variables. We assessed characteristics likely linked to environmental filtering, dispersal rates, drift, and species interactions. For environmental filtering effects, we used a subset of landform and climate related variables derived from the original datasets used by Bertin et al. 19 and Pfeiffer et al. 40 , and which can be found in Supplementary Table S5 online. Only non-redundant variables with apparent influence on the diversity parameters were considered. These were identified by performing partial least square (PLS) regressions on the different diversity metrics. Only environmental variables that contributed strongly to the retained PLS components (defined using a leave-one out cross-validation procedure and by identifying the first local minimum of the root-mean-squared error 75 ) were then retained for subsequent analyses. They included altitude, slope, mean annual wind speed (MAWS), mean annual number of snow days (SnowNDays), mean productivity (NDVI), and productivity variation (δNDVI). We used graph-based connectivity metrics as variables related to dispersal rates. These were quantified for each of the wetlands in this study using GRAPHAB 1.0 76 , by considering all other wetlands occurring within a 20 km radius from the focal wetland according to the wetland map generated by Bertin et al. 20 . The connectivity measures calculated were node degree (Dg, the number of patches close to the focal patch 77 ), closeness centrality (CCe, the mean distance from the focal patch to all other patches of its component 78 ), eccentricity (Ec, the maximum distance from the focal patch to another patch of its component 79 ), and betweenness centrality (BC, the sum of the shortest paths through the focal patch 76,80 ). We used wetland size and species co-occurrence as proxies for drift risk and plant species interactions, respectively. Wetland size was assessed based on remotely-sensed NDVI derived from LandSat 8 OLI satellite imagery (see Bertin et al. 20 ). Species co-occurrences were analyzed by calculating an average C-score for each wetland from the five plots. The C-score metric measures the level of segregation of all species pairs, which increases as average co-occurrences between all of the species pairs in the matrix decrease 81 . Thus, to allow an intuitive interpretation, the effects that we report for species co-occurrence are the additive opposite of the C-score effects. The C-score values were calculated with the R-package EcoSimR 81 .
Statistical analyses. Spatial structure analyses of the plant community parameters. We estimated spatial autocorrelation components of the community parameters using Moran's Eigenvector Maps 58 as implemented in the R packages spdep 82 and adespatial 83 . This framework allows producing spatial predictors, the MEMs, which model autocorrelation patterns at different spatial scales. They are obtained by eigen-decomposition of a spatial weighting matrix coding for connections between sites (i.e. wetlands in this instance), which produces orthogonal eigenvectors maximizing the spatial autocorrelation computed by Moran's I 84 . When associated with positive eigenvalues, the MEMs represent positive autocorrelation and vice versa 13 . The spatial weighting matrix that we used to generate the MEMs considered connections among all wetlands for which strength was inversely proportional to their Euclidean distance. To estimate S + (x) and S − (x) of community indices, we followed the procedure described in Biswas et al. 14,15 , which involved summing the squared correlations between each diversity parameter and the positive and negative MEMs weighted by their associated eigenvalue.
Influence of ecological factors on plant community parameters. We investigated the influence of each ecological factor (wetland connectivity and size, species co-occurrence, and environmental variability) on community indices by carrying out PLS-SEM 85 . This is a comprehensive technique based on the use of a variable covariance matrix that identifies relationships and causality between variables 86 with minimum requirements regarding measurement scales, sample sizes, and residual distributions 85 . The measurement model (i.e. the outer model, see Supplementary Figure S1 online) served to define latent variables representing the four ecological factors (wetland size, wetland connectivity, environmental variability, species co-occurrence) from the 12 measured ecological variables 87 . Environmental variability, wetland size, and connectivity were defined according to a formative mode and species co-occurrence using a reflective mode (see Supplementary Figure S1 online). Inspections of the loadings of the ecological variables on the latent variables indicated that the latent variables for wetland size and connectivity increased as wetland size and connectivity increased, respectively. The inner or structural model was constructed to estimate the effects of each of the four ecological factors on plant community characteristics (see Supplementary Figure S1 online). Eight separate PLS-SEM models were run to investigate these effects on each of the community metrics. Confidence intervals and significance of path coefficients were evaluated based on 10,000 bootstrapping iterations. These analyses were performed using the semPLS R-packages 85 .
Effects of assembly driving factors on positive and negative spatial autocorrelation components of community indices. To investigate whether the assembly driving factors influenced autocorrelation patterns of the community indices, we tested the relationship between the effects of the ecological factors (as estimated by their path coefficients) and the autocorrelation levels of the diversity parameters. To account for relatedness between the community parameters, we first analyzed the effects of their covariation on S + (x) and S − (x) levels. The pairwise Pearson's correlation matrix of the plant indices was transformed through principal coordinate decomposition (PCoA 16,88 ) using the ecodist R package 89 . Then, we analyzed the effects of the PCoA eigenvectors so-produced on S + (x) and S − (x) separately by carrying out generalized linear models using the R leaps package 90 . In these analyses, the autocorrelation levels of each community parameter were used as the dependent variables and the eigenvectors as the predictors. We ran the regsubset function, which performs an automatic selection procedure to look for the subset of eigenvectors that best explained the autocorrelation variation of the community indices. The best-fitting model, identified through a performance analysis performed with the AICcmodavg package 91 and using the Akaike´s Information Criterion corrected for small sample size (AICc), was then tested for sig- www.nature.com/scientificreports/ nificance. Since these analyses did not reveal any impact of the relatedness between the community indices on their autocorrelation levels (see Supplementary Table S2 online), influences of the plant community parameters on S + (x) and S − (x) were analyzed through multiple linear regressions. The dependent variable of these analyses were the autocorrelation levels of each community parameter and the predictors the PLS-SEM path coefficients measuring the direct effects of the ecological factors on each community characteristic. For these analyses, we also performed an exhaustive model selection to search for the subset of ecological factors that best explained the autocorrelation levels using, as above, the regsubset function and based on the AICc criteria. Normality of the residuals was verified using Shapiro tests.

Data availability
All the data used in the analyses are provided in Supplementary Tables S3, S4  www.nature.com/scientificreports/